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In this work we investigate the classical ferromagnetic XY-model in two di- 
mensions subject to a symmetry breaking field which impose a Z2 symmetry to 
the system. We used the broad histogram method combined with microcanoni- 
cal simulations and finite-size scaling analysis to estimate the critical temperature 
and critical exponents of this system with little computational effort. In addition, 
we present a general procedure which makes possible to use the broad histogram 
method for continuous systems, also when the macroscopic quantities needed by the 
method cannot be obtained analytically. Our results are robust under the choice of 
four different pseudo-random number generators. 

PACS: 02.70.Lq, 02.50.-r, 05.50.+q 

I. INTRODUCTION 

Consider a XY-model with Hamiltonian 

H = - cos(<9 4 ~ °j) ~ h pYl cos (ph) , (1) 

(ij) > 

where 9i denotes the angle between the spin <?j, in site i, and the x-axis, J is a coupling constant, 
the first summation extends over all pairs of nearest neighbors and the second one over all sites. 
The second term in (0) corresponds to the energy associated with a symmetry breaking field, char- 
acterized by a positive integer p and a constant h p . This model is important since it may describe 
many real magnetic crystals subject to symmetry-breaking crystalline fields fl. Experimental 
crystalline fields corresponding to p — 2, 4 and 6 are expected. For p = 1 the model corresponds 
to a XY-model with an external magnetic field. For p = 2 we have a system with Z2 symmetry 
(which we call XY2). This latter case is specially interesting because: (i) It was studied by Jose et. 
al. who were unable to ascertain the small /12 structure of the phase diagram 0] , (ii) The critical 
temperature T c for this system, neither for h% = 1> is not known exactly, (iii) Recently, Resende 
and Costa || have founded differences of 2% in the value of T c depending on the pseudo-random 
number generator employed, i.e., T c = 1.529(5) with RAN2 [| or T c = 1.491(6) with RANLUX 
in the fourth luxury level || , and (iv) the quantities needed to implement the recently introduced 
"Broad Histogram Method" (BHM) for this system are hard to compute. 

In this paper we have two main aims: (i) to show how to implement the BHM when there are 
no analytical expressions for its quantities and, (ii) to use the BHM, together with microcanonical 
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simulations, to calculate the critical temperature and the critical exponents for the XY2 model 

>>>2 li- 

In the next section we review the implementation of the BHM for continuous systems and we 
describe how to apply it to the XY2 model and, in extend, to any continuous model. In section 
(jlll|) some details on the sampling process employed are presented. Section (IV) shows the results 
of the finite-size scaling analysis. Finally, we summarize some comments and conclusion in section 



II. THE BROAD HISTOGRAM METHOD 



The broad histogram method (BHM) |3|,[7| is a novel analysis technique introduced in 1996 
by de Oliveira et. al. in order to calculate the spectral degeneracy g(E) of a given system. 
Instead of investigating the form of the histogram of visits, the new method computes g{E) from 
the microcanonical averages of a special macroscopic quantity N(X, AE) defined below. This 
new analysis method can be used in conjunction with any sampling strategy used previously 
(microcanonical, canonical, multicanonical, etc.). Since the basic quantities are microcanonical 
averages, the method does not depend on the form of histogram of visits and sample sets of different 
runs can be joined together just by cummulating them in the same histograms. Furthermore, each 
sample contributes with a number of several digits to these averages, a much richer information 
that just increasing by unity some entry in the histogram of visits. As a result, the BHM can 
estimate thermodynamic averages with more precision on broader temperature ranges than other 
reweighting methods fi^|. 

Consider some symmetric relationship R connecting pairs of states in the phase space of the 
system and define a new macroscopic quantity iV(X, AE) as the number of connections starting 
from state X that changes the energy of this state by AE, Now, consider the subspaces T(E) 
and T(E + AE) of all states with energy E and E + AE, respectively. Since R is symmetric, the 
number of connections going from T(E) to T(E + AE) equals the number of connections coming 
back from T(E + AE) to T(E), or 

g(E)(N{E, AE)) = g{E + AE){N(E + AE, -AE)) , (2) 

where g(E) is the number of microstates with energy E and (N(E,AE)) is the average of 
N(X.,AE) on T(E). This is the so called broad histogram equation. It means that the ratio 
g(E + AE)/g(E) can be obtained by estimation of the averages (N(E,AE)). Equation (0) is 
completely general and is exact for any statistical model JnJ. Starting from an initial value g{Eo) 
and applying the method iteratively, the density of states function g(E) can be obtained in steps 
of AE on a whole energy range. 

By taking logarithm at both sides of eq. (|J) we arrive to the following expression for 
P(E):=d\n[g(E)}/dE: 

m * ^ In [ {N{E + AEj _ AE)) ) ■ (3) 

Then, an alternative procedure is to use this expression to estimate @(E) at several energy values 
and integrate it to obtain an approximation of S{E):=hx[g{E)], the microcanonical entropy. 

The BHM is an analysis procedure on a set of samples obtained previously. The samples itself can 
be colleted with every sa mpling method avaible as, for instance, microcanonic al j ll[ Jl2]|, canonical 
T3|.[l4f (including cluster p5 16 algorithms) or multicanonical [p^| - [22[ and 1/fc [^3| procedures, and 



some of them have been already tested [[| p[|2j| |3l[ . The averages {N(E,AE)) can be estimated 
just by cumulating the quantities iV (X, AE) on histograms in E and by dividing them into the 
number of samples at each energy bin, so far states with the same energy are taken into the sample 
with the same probability (This is the usual procedure to compute microcanonical averages). 

If we are interested in canonical averages (Q)t a t temperature T of a certain quantity Q, we 
need an additional histogram in E to cumulate Q{X) for all states in the sample and compute 
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the microcanonical averages (Q(E)). Once g(E) and (Q(E)) have been determined, the canonical 
aveage {Q)t of a certain quantity Q are computed as 

(QW = ^Y.iQiE)) g(E)eM~E/k B T] 

Z(T) = J29(E)exp[-E/k B T] (4) 

E 

for any desired temperature, where fee is the Boltzmann constant. 

The BHM method was originally used to investigate discrete models but .recently, Muhoz and 
Herrmannn have shown how to extend it for continuous systems, as the XY [^,^9) and Heisenberg 
p8| models. The relationship R is replaced by a symmetric transition probability i?(Y|X) = 
i?(X|Y), and the number of connections -ZV(X, AE) is redefined as the probability, /ab(X, AE), 
to obtain an energy change between AE and AE + dAE with this transition probability R. Hence, 
the problem reduces to that of calculating the probability density function (p.d.f.) /a_e(X, AE) 
for a given configuration X. 

This calculation can be easily done for the XY and Heisenberg models, where a partial analytical 
computation is easy to perform. For the XY model with Z2 symmetry the problem is more difficult. 
Consider the transition probability R defined by 

{T^rjy- if Y can be reach from X by a 
single-spin turn . (5) 

otherwise 

Here, N is the number of sites. This transition probability corresponds to the following protocol 
of virtual random movements: first, a site i is chosen at random (with probability 1/iV); then, a 
new angle 9i new is chosen at random on the interval [— 7r, tt). Thus, 9i new is a random variable 
with p.d.f. 

c ia\ 1 V^ 71 " : if vr < 9 t now < 7r , . 

M„wW = | o : otherwise ' (6) 

If the spin i was chosen to be turn, its bonds and interaction with the external field are the only 
ones that change the energy value. Let us define the energy associated to spin i as 

7 = - £ cos (^ cos ( 2 ^) 
i 

= —Ai cos(0j) — Bi sin(#j) — — cos(2#i) 

^giA^BiA), (7) 

where Aj=Y]j cos(#j), Bj=Y]j sm(#j) and the summations run over all nearest neighbors of spin i. 
If 6i changes to 9i now the energy ej changes to ej now . What is the p.d.f. of this random variable? 
To find the p.d.f. / £ ™w(e) for a specific e we need to solve the equation 

g(Ai,Bi,e)- e = 0. (8) 
Denoting its real roots by 9i,02, the p.d.f. is written as p3] 

J ^ y> \dg/d9(A i ,B i ,0 1 )\^ \dgldO(Ai,BMC \dg/d9{A ll B l ,9 n )[ {) 

The difficulty here is to find the roots $i , 6%, ■ ■ ■ , 6 n . The equation (|^) is fourth order in 9 (con- 
sequently the maximum value for n is 4) and there are analytical expressions to compute the 
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roots, but they are too CPU-time demanding to be computed for every site i. Moreover, the exact 
calculation of these roots for other models could be impossible. 

For these reasons we propose to calculate the roots numerically and to store a table 
f eiMew (Ai,Bi,e) for some values of Ai, Bi and e. Other values not stored in the table can be 
founded by interpolation techniques. This approach is very fast, since we need only to calculate 
Ai, Bi and E\ for each selected spin and to read from the table the value for the respective p.d.f.. 
Figure (1) shows the p.d.f. f £i (£i) for some values of Ai and Bi, obtained from a table with 128 
entries for Ai, 128 for Bi and 256 for e. Despite of these few values we did not find any significant 
difference between using linear interpolation or taking the nearest value in the table. Furthermore, 
nearly the same results are obtained with only 128 entries for £j. 

Because AE=Ei ncw —ei, its p.d.f. for a change in site i is f^ E (AE)=f eiaew (Ai, B,,e, + AE). 
Finally, since all sites are equally probable, for a given configuration X we obtain 

f AE (X, AE) = ^J2 A. .™ (A, B u e l + AE). (10) 

i 

where the sum is over all sites. 



III. SAMPLING PROCEDURE 



Toward a finite-size scaling analysis of this two-dimensional model we investigated system sizes 
of L=6, 8, 10, 16, 20, 40, 60 and 100 with periodic boundary conditions. A single value of A_E=4.0 was 
chosen for the BHM analysis. Since we are interested only in thermodynamic functions near the 
critical temperature it is not necessary to sample the whole energy range of positive temperatures. 
Instead, the following schema was adopted: For each system size we set a maximal and a minimal 
energy values -E max an d -E^in and chose NE equally-spaced energy points between them. Our 
goal is to estimate /3(E) at each point E, as 

HK ' AE b \(N(E + AE b /2,~AE))J y bJ ' v ' 

and to integrate it by the trapezoidal rule. 

Samples were taken by performing microcanonical simulations on small windows of width AEb, 
that choose and turn one spin at random and accept the new state only if it falls into the window. 
This procedure takes all states inside the window with the same probability. Then, most of the 
samples will cumulate on the high-energy border of the window, because the density of states is a 
fast-growing function. Therefore, (N(E + AEb/2, —AE)) is well estimated from a microcanonical 
simulation in a window of width AEb centered at E, and (N(E — AEb/2, AE)) from a micro- 
canonical simulation in a window of width AEb centered at E — AEb. A schematic description of 
this sampling procedure is shown in figure (||) . Table | shows the values of E max E m i n and NE 
used for each system size. 

Summarizing, we ran two microcanonical simulations for each value of (3(E). The microcanonical 
averages (M)e, (M 2 ) e and (M 4 ) e were computed at each energy point by averaging the M- values 
in both windows. Inside each window 1000 MCSS were discarded before sampling and, then, 100 
samples spaced by 10 MCSS were taken. The whole procedure was repeated 20 times and the 
data of all 20 runs were cumulated in the same histograms. 



IV. FINITE-SIZE SCALING ANALYSIS 

We performed a finite-size scaling analysis J33|,[34| to estimate the critical properties of this model. 
We have divided our analysis in two parts; first, the determination of the critical temperature and, 
second, the estimation of the critical parameters. 
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A. Critical temperature and Binder cumulant 



One of the questions we want to elucidate is the dependence of the results on the random 
number generator employed. Hence, we repeated our simulations with four different random 
number generators, i.e., RAN2 §j, RANLUX with the fourth luxury level [§, RANIBM [j§ and a 
lagged-FIBONACCI ]35| , |36| ] random number generator with primitive trinomial (110503,53719,0) 
jjTJ, initialized with RAN2. 

Figure (^) shows the fourth-order Binder cumulant 

Ul{t) - 1 -^w% (12) 

computed for all system sizes with the four random number generators. By looking at the limits of 
the smallest rhombus crossed by all curves we define maximal and minimal values for the critical 
coupling (3t c and the universal parameter U* := Ul(T c ) (table |ll|). The most precise results are 
obtained with the FIBONACCI and the less ones with RANLUX. By combining the four results 
through a maximal likehood estimator |§ we obtain T c = 1.5329(4) (j3 c = l/k B T c = 0.6524(2)) 
and U* = 0.6118(6) , in good agreement with the value U* = 0.618 for the 2D-Ising universality 
class f39f (here and in the following the number between parenthesis represents a 2a error bar, 
with a confidence level of 95%. For instance, 2.158(24) = 2.158 ± 0.024 and a = 0.012). 

To corroborate this value for the critical temperature we have plotted the temperatures of 
maximal heat capacity Tc v (L) and maximal magnetic susceptibility Tx s {L) as a function of the 
system size L. In the thermodynamic limit 

T Cv {L)=T c + aL- 1/v ; T Xs (L) = T c + bL- 1/v . (13) 

Since v = 1.0 for the 2D-Ising universality class, we expect a linear behavior of Tc v {L) and Tx s {L) 
against 1/L. Figure (^) shows this linear behavior with the four random number generators. The 
linear regressions were obtained from the data points at L = 10, 16, 20, 40, 60 and 100. All figures 
coincide to estimate a critical temperature near T c — 1.53, in excellent agreement with our estimate 
from the Binder cumulant and in concordance with the value of T c = 1.529(5) obtained in ref. || 
with RAN2, but far away of value T c = 1.491(6) obtained in the same reference with RANLUX. 



B. Critical exponents 

To compute the critical exponents of this model we used the data obtained with the FIBONACCI 
random number generator, that has showed to be the most precise one. 

The critical exponent v can be obtained by two different procedures. First, it can be estimated 
from the previous data for the Binder cumulant by using the relation |34| 



dU L 



dT 



-L~ x l v . (14) 



Figure (|a) shows the power-law fit for this relation at T c = 1.5329. The value l/v = 1.009(18) 
(y = 0.991(16)) obtained is consistent with the value v — 1.0 of the 2d-Ising universality class. 
Second, it can be estimated from ( |l3| ) by plotting Tx s (L) — T c against 1/L (figure (^b)). This 
second value of l/v = 0.920(24) (v = 1.087(28)) is not so precise as the first one. By combining 
both results in a maximal likehood estimator, we obtain v — 1.015(16). 

Figure (^a) shows the average magnetization per site {M)t c at T c = 1.5329 as a function of the 
system size L, in order to compute the critical parameter (3 /v from the finite-size behavior 

±(M) Te ~L-f > /». (15) 
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The value of j3/v = 0.1340(34) is higher than the value 1.25 of the 2d-Ising universality class. 
This can be a small-sizes effect, since the same estimation excluding L = 6,8, 10 and 16 gives 
j3Jv = 0.1285(20). By excluding in consecutively the system sizes L = 6, 8, 10 and 16 we obtained 
different values of fi/v. These values, when plotted against X/L, show a linear behavior (L is 
the minimal system size taken into account) (figure (^]b)). The linear extrapolation gives a final 
estimation for (3/v = 0.126(2). Figure (||c) shows the FSS collapse of the average magnetization, 
by taking v = 1.015 and j3/v = 0.126. 

To compute j/u we have used the three finite-size scaling relations 



±Xs(T c )~L^ , (16) 

^%ax~L^ , (17) 

-l(Af 2 ) Tc ^L^ . (18) 



Figure (Qa) shows the results form these three methods. The results are 7/1/ = 1.745(8) (eq. 17), 
7/1/ = 1.741(4) (eq. [t8|) and j/u = 1.798(4) (eq. By excluding consecutively the system sizes 
L = 6, 8, 10 and 16 we obtain the values shown in figure (0b). The lar ger variation corresponds 
to the estimations from Xs max , as expected The linear extrapolation of these values does 

not coincide for L — > 00, Nevertheless, they all intersect around 7/1/ = 1.742(4), when the system 
sizes L = 6, 8 and 10 have been excluded. Furthermore, it coincides with the maximal-likehood 
estimator of the three values, 7/1/ = 1.742(4), in good agreement with the value for the 2d-Ising 
universality class, 7/1/ = 1.75. Figure (Qc) shows the FSS collapse of the magnetic susceptibility 
with v = 1.015 and j/u = 1.742. 

The more difficult parameter to be directly computed is a, the critical parameter for the diver- 
gence of the specific heat. On one side, assuming hyperscaling a = 2 — du to be valid and setting 
our estimate of v = 1.015(16) we obtain for d — 2 (the lattice dimension) a = —0.030(32), which 
corresponds to a finite cusp, 

C w (£) = C W (L = 00) - aL a/v . (19) 

On the other side, for the 2d-Ising universality class a — 0.0 and Cv max shows a logarithmic 



divergence 45 



Cv max (L) = alnL + b + 0[lnL/L] . (20) 
Furthermore, one predicts the scale behavior p6| , ^3|] 

C VL (x) = A lnL + Q^(x)+O[lnL/L] , (21) 
where x ~ L(T — T c )/T c . The constant Aq can be deduced from the scaling behavior of C VL (T C ) 

Cv L (T c ) = A hxL + Q( '(0) + 0[ln£/i] . (22) 

Figure (Ha) shows both a power-law fit with a/v = -0.039 (C„ majE (L) = 25.0 + 3.200(4) x 
L -o.03938(iV)) and a i g ar i t hmic fit (C Vmiuc (L) = 0.62(10) + 0.853(31) InL) for C Vm ^(L) (circles). 
The value of C Vmta .(oo) — 25.0 for the first one was set by looking for the larger correlation 
coefficient r of the power-law fit. This maximal value r = 0.9989758 changes, however, only 
in the fifth decimal place for 16.0 < CV, max (oo) < 297.0, corresponding to —0.0732(28) < a < 
0.00281(10). We take, hence, a = —0.04(8), in good agreement with the value from hyperscaling. 
The logarithmic fit, with correlation coefficient r — 0.99888879 is almost indistinguishable from 
the power-law fit. 

Figure (ga) also shows a logarithmic fit for C VL (T c ) (squares). This fit estimates Aq = 0.818(20) 
from (p2) with correlation coefficient r = 0.9995723. Figure (&>) shows the FSS collapse of the 
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specific heat according with (£L|) for this value of A$. All system sizes collapse well, except for 
L = 100 at T < T C) due, perhaps, to a lack of simulational points for this system size. All these 
results support the thesis of a logarithmic anomaly of the specific heat for this model. 

Our results are summarized in table They are in good agreement with the expected values 
for the 2d-Ising universality class and [j3j. 



V. COMMENTS AND CONCLUSIONS 



We have constructed a BHM analysis procedure for the XY-model with Zi symmetry. With 
this example, we have shown how to construct a BHM analysis when the quantities N(E, AE) 
cannot be easy computed. Our strategy is to construct a table for the function /* (eq. [)]), where 
the indexes are given by the variables which define the energy per spin £j, and the values of the 
function are determined by numerical methods. Other values not stored in the table are found by 
interpolation techniques. This approach is very fast, since only a memory access, an interpolation 
and a sum are needed to compute N(E,AE). It is only restricted by the amount of memory 
available. This approach is general and can be extended to other complex models. 

We have employed the BHM analysis procedure together with microcanonical simulations to 
perform a finite-size scaling analysis of the XY-model with Zi symmetry, for system sizes 6 < L < 
100. Our calculations give a very precise estimates for the critical temperature, T c = 1.5329(4) 
and of the Binder cumulant U* = 0.6118(6). Furthermore, our results for T c and U* are almost the 
same for the four pseudo-random number generators employed, RAN2, RANLUX, FIBONACCI 
and RANIBM, with the most precise results for FIBONACCI and the least ones for RANLUX. 
Our estimations for the critical exponents arc in good agreement within the error bars with the 
theoretical predictions of the 2d-Ising universality class and with other Monte Carlo studies. 

Our value of T c agrees with the result T c = 1.529(5) from ||| with RAN2 and clearly disagrees 
with the result T c = 1.491(6) from the same reference with RANLUX. To clarify this point, 
figure (||) shows our results for the magnetic susceptibility in a 20 x 20 system, obtained with 
different random number generators and including also some reference points calculated with 
the Metropolis algorithm, which reproduce the figure (1) of ref. H. In that previous work, the 
Metropolis algorithm with RANLUX showed a peak at a clearly lower temperature (T = 1.55) 
than the same algorithm with RAN2 (T = 1.58). In contrast, all our data coincide to show a peak 
around T — 1.58, including the Metropolis estimations with RANLUX. In our opinion, the result 
of T c = 1.491(6) could be due to an incorrect implementation of RANLUX. 

The computational time of the BHM analysis increases proportional to N, the number of sites, 
in the same proportion as the time of one MCSS with the microcanonical sampling procedure. A 
special care should be taken, however, to assure that enough microcanonical simulations contribute 
to each temperature average. According to eq. (^) , the contribution of a sample with energy E to 
the average {Q)t is weighted by the factor g{E) exp(—E(x)/kBT), which distributes approximately 
as a Gaussian with center at (E)t and width a e(T)=T 'yj 'k-g,C v {T) . Therefore, all samples in the 
range (E)t — 3<je < E < (E)t + 3<te will contribute to (Q)t- For small system sizes, the whole 
range of energies contributes. For large ones, we have tried to maintain a constant number of 
simulations inside 1o\e. The number of simulations per energy interval can be chosen by assuming 
C v (T)~L d for a < and C V (T) ~ i 1 d + a /' J otherwise. Table IV shows the number of microcanonical 
simulations in lcr B at x = L(T — T c )/T c = —1 by assuming C v (x = —1) = 0.8181nL — 0.5 (see 
figure (§>)). 

In conclusion, we have shown that the BHM method can be used to obtain quality results from 
the finite-size scaling analysis of a non-trivial system, also if the quantities required to implement 
the method are hard to compute analitically. These results support the usefulness and generality 
of the BHM. 
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L 


6 


8 


10 


16 


20 


40 


60 


100 


^max 


0.0 


0.0 


0.0 


0.0 


-0.35 


-0.35 


-0.35 


-0.35 


J-^min 


-1.0 


-1.0 


-1.0 


-1.0 


-0.70 


-0.70 


-0.70 


-0.70 


NE 


200 


200 


200 


320 


140 


280 


420 


700 



TABLE I. Parameters of the microcanonical simulations. 



Generator 


RAN2 


RANLUX 


FIBONACCI 


RANIBM 


Combined 




1.5317(16) 


1.5306(27) 


1.5330(4) 


1.5328(13) 


1.5329(4) 


U* 


0.6110(29) 


0.6139(32) 


0.6118(6) 


0.6109(42) 


0.6118(6) 



TABLE II. Critical temperature T c and Binder cumulant U* for the two-dimensional XY-model with 
Z2 symmetry. 





Generator 


V 


0/v 


7/" 


a/v 


This study 


FIBONACCI 


1.015(16) 


0.126(2) 


1.742(8) 


-0.04(8) 


Ref. [ 


! 


RANLUX 




0.14(6) 


1.71(7) 


0.00(3) 


Ref. 


i 


RAN 2 




0.13(4) 


1.75(9) 


0.05(4) 


2d Ising 




1.000(0) 


0.125(0) 


1.750(0) 


0.00(0) 



TABLE III. Critical Parameters for the XY-model with Z2 symmetry. 



L 


6 


8 


10 


16 


20 


40 


60 


100 


NE/a E 


14 


12 


11 


13 


14 


16 


17 


18 



TABLE IV. Number of microcanonical simulations at x = — 1 for the settings of Table 
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FIG. 1. Probability distribution for : A = 0, B = 4 (dot-dashed line), A = 2.25, B = 2.25 (continuous 
line), A = -1.75, S = 2.25 (dashed line) and A = 0.5625, B = -0.1875 (long-dashed line) 



FIG. 2. Sampling procedure. 



FIG. 3. Binder cumulant for the XY-model with Z2 symmetry obtained with (a) the RAN2, (b) the 
RANLUX, (c) the FIBONACCI and (d) the RANIBM random number generators. 

FIG. 4. Size dependence of the location of the peak for the magnetic susceptibility (circles) and the 
heat capacity (squares) obtained with (a) the RAN 2, (b) the RANLUX, (c) the FIBONACCI and (d) 
the RANIBM random number generators. 

FIG. 5. Estimations of u from the slope of the Binder cumulant Ul at T c (a) and the position of the 
peak in the magnetic susceptibility T Xs (b). We adopted T c = 1.5329. 

FIG. 6. (a) Estimations of j3/v from the average magnetization per site at T c — 1.5329. (b) Linear 
extrapolation of this estimations, (c) FSS of the average magnetization with v — 1.015 and fi/v — 0.126 
for system sizes L = 16, 20, 40, 60 and 100. 

FIG. 7. (a) Estimations of 7/1/ from the magnetic susceptibility per site at T c — 1.5329 (circles), the 
maximal value of magnetic susceptibility (squares) and the average value of the square of the magnetization 
at T c (diamonds), (b) Linear extrapolation of this estimations, (c) FSS of the magnetic susceptibility 
with v = 1.015 and 7 = 1.742 for system sizes L = 16, 20, 40, 60 and 100. 

FIG. 8. (a) Power-law fit (dashed line) and logarithmic fit (solid line) of the maximal specific heat 
(circles) as a function of the system size L, and logarithmic fit (dot-dashed line) of the specific heat at 
T c = 1.5329 (squares), (b) Finite size scaling of the specific heat assuming a logarithmic anomaly for 
system sizes L = 16, 20, 40, 60 and 100. 

FIG. 9. Magnetic susceptibility for a 20x20 system, computed from BHM on pairs of microcanonical 
simulations (BHM-IN) and from Metropolis simulations (MC), both using several pseudo-random number 
generators. 
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